#
# Plot the histogram produced by genomeCoverageBed as from
# $ genomeCoverageBed -max 200 -ibam %(bam)s -g %(bedtools_genome)s > %(bam)s.genomeCoverageBed.hist
#
# Note: All outputs in dir bam_clean ending with .bam.genomeCoverageBed.hist are plotted
#
# USAGE: R CMD BATCH ~/python-packages/chipseq-pipeline/plot_pileup_histograms.R

hist_table<- list.files('bam_clean')
hist_table<- grep('\\.hist$', hist_table, perl= TRUE, value= TRUE)

w= 900
h= (w / 4) * ceiling(length(hist_table)/4)
png('redmine_wiki/pileup_histograms.png', width= w, height= h, pointsize= 18)

par(mfrow= c(ceiling(length(hist_table)/4), 4), mar= c(4,4,1.5,0.1), las= 1, mgp= c(2, 0.75, 0), oma= c(0,0,2,0))
for(hfile in hist_table){
    name<- sub('.bam.genomeCoverageBed.hist', '', hfile)
    ghist<- read.table(paste('bam_clean', hfile, sep= '/'), sep= '\t', stringsAsFactors= FALSE)
    ghist_genome<- ghist[ghist$V1 == 'genome' & ghist$V2 > 0,]
    y= log2(ghist_genome$V3)
    plot(y= y, x= ghist_genome$V2, type= 'n', xlab= 'No. reads (pileup)', ylab= 'Log2 No. of bases', main= name, cex.main= 0.85, ylim= c(0, max(y)))
    grid(col= 'lightblue')
    points(y= y, x= ghist_genome$V2, type= 'l', lwd= 1.5, col= 'firebrick4')
    }
mtext(text= 'Histograms of depth of sequencing', outer= TRUE, side= 3, line= 0.75, font= 2)
dev.off()